Libraries required
library(limma)
library(edgeR)
library(plgINS)
library(sva)
library(ggplot2)
library(dplyr)
library(patchwork)
library(SummarizedExperiment)
library(reshape2)
library(autoplotly)
library(pheatmap)
library(viridis)
library(RColorBrewer)
library(ggplotify)
Load salmon objects
load("input/SC_controls_rnaseq_salmon.tds.RData")
lab <- salmon
load("input/SC_controls_rnaseq_lit_salmon.tds.RData")
lit <- salmon
counts_lab <- lab@tx.counts
counts_lit <- lit@tx.counts
counts_lab <- counts_lab[,grep(pattern = "Adult", x = colnames(counts_lab))]
counts_lab <- data.frame(Genes = rownames(counts_lab), counts_lab)
counts_lit <- counts_lit[,grep(pattern = "PNW8", x = colnames(counts_lit)), drop = FALSE]
counts_lit <- data.frame(Genes = rownames(counts_lit), counts_lit)
counts <- merge(counts_lab, counts_lit)
rownames(counts) <- counts$Genes
counts <- counts[,-1]
pdata <- data.frame(Samples = colnames(counts), Group = c(rep("PNW20", 6), "PNW8"))
Differeitial analysis using limma
PNW8 vs PNW20
design <- model.matrix(~ 0 + Group, data = pdata)
colnames(design) <- gsub(pattern = "Group", replacement = "", x = colnames(design))
y <- DGEList(counts = counts)
keep <- filterByExpr(y, design, min.count = 10)
y <- y[keep, ]
dds <- calcNormFactors(y)
v <- voom(dds, design = design)
contrast.matrix <- makeContrasts(PNW20 - PNW8, levels = design)
fit <- lmFit(v)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
pnw8.pnw20 <- as.data.frame(topTable(fit2, coef = 1, number = Inf))
pnw8.pnw20 <- data.frame(Genes = rownames(pnw8.pnw20), pnw8.pnw20, stringsAsFactors = F)
pnw8.pnw20_sig <- pnw8.pnw20[abs(pnw8.pnw20$logFC) >= 1 & pnw8.pnw20$adj.P.Val <= 0.05, ]
cpm Counts data and pData
y <- DGEList(counts = counts)
dds <- calcNormFactors(y)
cpm <- cpm(dds, log = T)
data <- list(cpm = cpm, pData = pdata)
save(data,
file = "./output/data_pData.RData", compress = T,
compression_level = 3
)
Results
dea.list <- list(
`pnw20 vs pnw8` = as.DEA(pnw8.pnw20)
)
dea.limma <- list(
`pnw20 vs pnw8` = pnw8.pnw20
)
Save RData files
save(dea.list,
file = "./output/dea_SC_Controls_pnw8_pnw20.DEA.RData", compress = T,
compression_level = 3
)
save(dea.limma,
file = "./output/limma_SC_Controls_pnw8_pnw20.RData", compress = T,
compression_level = 3
)
writexl::write_xlsx(x = dea.limma, path = "output/dea_results_pnw8_pnw20.xlsx", col_names = T, format_headers = T)
MA plots
res <- pnw8.pnw20
res$threshold <- as.factor(res$adj.P.Val < 0.05)
p1 <- ggplot(data = res, aes(
x = res$AveExpr,
y = res$logFC,
colour = threshold
)) +
geom_point(alpha = 0.5, size = 1.8) +
geom_hline(aes(yintercept = 0), colour = "blue", size = 1) +
ylim(c(
-ceiling(max(abs(res$logFC))),
ceiling(max(abs(res$logFC)))
)) +
ggtitle("PNW20 vs PNW8") +
labs(subtitle = "loess fit") +
xlab("Mean expression") +
ylab("Log2 Fold Change") +
theme(
plot.title = element_text(face = "bold", size = 20, hjust = 0.5),
plot.subtitle = element_text(size = 16, hjust = 0.5, face = "italic", color = "blue"),
axis.title.x = element_text(face = "bold", size = 15),
axis.text.x = element_text(face = "bold", size = 12),
legend.title = element_text(face = "bold", size = 15),
legend.text = element_text(size = 14)
) +
scale_colour_discrete(name = "p.adjusted < 0.05") +
stat_smooth(se = FALSE, method = "loess", color = "red", formula = y ~ x, size = 1)
## Warning: Use of `res$AveExpr` is discouraged. Use `AveExpr` instead.
## Warning: Use of `res$logFC` is discouraged. Use `logFC` instead.
## Warning: Use of `res$AveExpr` is discouraged. Use `AveExpr` instead.
## Warning: Use of `res$logFC` is discouraged. Use `logFC` instead.
Heatmap of most variable genes
xvar <- apply((data$cpm + 1), 1, var)
genes500 <- head(sort(xvar, decreasing = TRUE), n = 500)
x500 <- data$cpm[rownames(data$cpm) %in% names(genes500), ]
pheatmap(mat = x500, scale = "row", show_rownames = F, col = viridis(100),
show_colnames = T,
main = "Clustering of samples based on top 500 variable genes")
PCA
pca <- prcomp(t(data$cpm))
# ggplotly(
autoplot(pca,
data = data$pData, colour = "Group",
frame = TRUE, frame.type = "norm", size = 2
) +
ggtitle("") +
theme(axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20),
plot.title = element_text(size = 25, face = "bold"),
legend.text=element_text(size=20),
legend.title=element_text(size=20, face = "bold"))
## Too few points to calculate an ellipse
Volcano plots
tab <- lab@tx.info
tab <- tab[,c(4,6,8,12)]
rownames(tab) <- NULL
tab <- tab[!duplicated(tab),]
tab <- merge(tab, pnw8.pnw20, by.x = "transcript", "Genes")
tab$Significant <- "Not significant"
tab$Significant[tab$adj.P.Val <= 0.05 & abs(tab$logFC) >= 1] <- "PNW8"
tab$Significant[tab$Significant == "PNW8"] <- ifelse(tab$logFC[tab$Significant == "PNW8"] > 0, "PNW20", "PNW8")
tab$Significant[tab$Significant == "PNW8"] <- "Ribo0"
tab$Significant[tab$Significant == "PNW20"] <- "polyA"
a <- ggplot(tab, aes(x = logFC, y = -log10(adj.P.Val))) +
geom_point(aes(color = Significant), size = 2, alpha = 0.5) +
scale_color_manual(values = c("black", "#0072B5FF", "#BC3C29FF")) +
xlim(
-(ceiling(max(abs(tab$logFC)))),
ceiling(max(abs(tab$logFC)))
) +
ylim(0, max(-log10(tab$adj.P.Val), na.rm = TRUE) + 1) +
xlab(bquote(~ Log[2] ~ "FC")) +
ylab(bquote(~ -Log[10] ~ adjusted ~ italic(P))) +
labs(
title = "PNW20 vs PNW8"
) +
theme_bw(base_family = "Helvetica") +
geom_vline(
xintercept = c(-1, 1),
linetype = "longdash",
colour = "black",
size = 0.4
) +
geom_hline(
yintercept = -log10(0.05),
linetype = "longdash",
colour = "black",
size = 0.4
)
a1 <- a + facet_wrap(~type, scales = "free_y")
a2 <- a + facet_wrap(~type2, scales = "free_y")
References
report::cite_packages(session = sessionInfo())
## Warning in citation(pkg_name): no date field in DESCRIPTION file of package
## 'plgINS'
SessionInfo
devtools::session_info() %>%
details::details()
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 3.6.2 (2019-12-12)
os Ubuntu 16.04.7 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Europe/Zurich
date 2020-12-11
─ Packages ───────────────────────────────────────────────────────────────────
package * version date lib
annotate 1.64.0 2019-10-29 [1]
AnnotationDbi 1.48.0 2019-10-29 [1]
assertthat 0.2.1 2019-03-21 [1]
autoplotly * 0.1.2 2018-04-21 [1]
backports 1.1.8 2020-06-17 [1]
base64enc 0.1-3 2015-07-28 [1]
Biobase * 2.46.0 2019-10-29 [1]
BiocGenerics * 0.32.0 2019-10-29 [1]
BiocManager 1.30.10 2019-11-16 [1]
BiocParallel * 1.20.1 2019-12-21 [1]
Biostrings 2.54.0 2019-10-29 [1]
bit 4.0.4 2020-08-04 [1]
bit64 4.0.2 2020-07-30 [1]
bitops 1.0-6 2013-08-17 [1]
blob 1.2.1 2020-01-20 [1]
bookdown 0.20 2020-06-23 [1]
callr 3.4.3 2020-03-28 [1]
checkmate 2.0.0 2020-02-06 [1]
cli 2.0.2 2020-02-28 [1]
cluster 2.1.0 2019-06-19 [1]
colorspace 1.4-1 2019-03-18 [1]
crayon 1.3.4 2017-09-16 [1]
data.table 1.13.0 2020-07-24 [1]
DBI 1.1.0 2019-12-15 [1]
DelayedArray * 0.12.3 2020-04-09 [1]
desc 1.2.0 2018-05-01 [1]
DESeq2 1.26.0 2019-10-29 [1]
devtools 2.3.1 2020-07-21 [1]
digest 0.6.25 2020-02-23 [1]
dplyr * 1.0.1 2020-07-31 [1]
edgeR * 3.28.1 2020-02-26 [1]
ellipsis 0.3.1 2020-05-15 [1]
evaluate 0.14 2019-05-28 [1]
fansi 0.4.1 2020-01-08 [1]
farver 2.0.3 2020-01-16 [1]
foreign 0.8-76 2020-03-03 [1]
Formula 1.2-3 2018-05-03 [1]
fs 1.5.0 2020-07-31 [1]
genefilter * 1.68.0 2019-10-29 [1]
geneplotter 1.64.0 2019-10-29 [1]
generics 0.0.2 2018-11-29 [1]
GenomeInfoDb * 1.22.1 2020-03-27 [1]
GenomeInfoDbData 1.2.2 2019-11-18 [1]
GenomicRanges * 1.38.0 2019-10-29 [1]
GEOquery 2.54.1 2019-11-18 [1]
ggfortify 0.4.10 2020-07-27 [1]
ggplot2 * 3.3.2 2020-06-19 [1]
ggplotify * 0.0.5 2020-03-12 [1]
glue 1.4.1 2020-05-13 [1]
gridExtra 2.3 2017-09-09 [1]
gridGraphics 0.5-0 2020-02-25 [1]
gtable 0.3.0 2019-03-25 [1]
Hmisc 4.4-1 2020-08-10 [1]
hms 0.5.3 2020-01-08 [1]
htmlTable 2.0.1 2020-07-05 [1]
htmltools 0.5.0 2020-06-16 [1]
htmlwidgets 1.5.1 2019-10-08 [1]
httr 1.4.2 2020-07-20 [1]
insight 0.9.0 2020-07-20 [1]
IRanges * 2.20.2 2020-01-13 [1]
jpeg 0.1-8.1 2019-10-24 [1]
jsonlite 1.7.0 2020-06-25 [1]
knitr 1.29 2020-06-23 [1]
labeling 0.3 2014-08-23 [1]
lattice 0.20-41 2020-04-02 [1]
latticeExtra 0.6-29 2019-12-19 [1]
lazyeval 0.2.2 2019-03-15 [1]
lifecycle 0.2.0 2020-03-06 [1]
limma * 3.42.2 2020-02-03 [1]
locfit 1.5-9.4 2020-03-25 [1]
magrittr 1.5 2014-11-22 [1]
Matrix 1.2-18 2019-11-27 [1]
matrixStats * 0.56.0 2020-03-13 [1]
memoise 1.1.0.9000 2020-05-06 [1]
mgcv * 1.8-31 2019-11-09 [1]
munsell 0.5.0 2018-06-12 [1]
nlme * 3.1-148 2020-05-24 [1]
nnet 7.3-14 2020-04-26 [1]
patchwork * 1.0.1 2020-06-22 [1]
pheatmap * 1.0.12 2019-01-04 [1]
pillar 1.4.6 2020-07-10 [1]
pkgbuild 1.1.0 2020-07-13 [1]
pkgconfig 2.0.3 2019-09-22 [1]
pkgload 1.1.0 2020-05-29 [1]
plgINS * 0.1.5 2020-07-10 [1]
plotly 4.9.2.1 2020-04-04 [1]
plyr 1.8.6 2020-03-03 [1]
png 0.1-7 2013-12-03 [1]
prettyunits 1.1.1 2020-01-24 [1]
processx 3.4.3 2020-07-05 [1]
ps 1.3.3 2020-05-08 [1]
purrr 0.3.4 2020-04-17 [1]
R6 2.4.1 2019-11-12 [1]
RColorBrewer * 1.1-2 2014-12-07 [1]
Rcpp 1.0.5 2020-07-06 [1]
RCurl 1.98-1.2 2020-04-18 [1]
readr 1.3.1 2018-12-21 [1]
remotes 2.2.0 2020-07-21 [1]
report 0.1.0 2020-03-19 [1]
reshape2 * 1.4.4 2020-04-09 [1]
rlang 0.4.7 2020-07-09 [1]
rmarkdown 2.3 2020-06-18 [1]
rmdformats 0.4.0 2020-06-07 [1]
rpart 4.1-15 2019-04-12 [1]
rprojroot 1.3-2 2018-01-03 [1]
RSQLite 2.1.4 2019-12-04 [1]
rstudioapi 0.11 2020-02-07 [1]
rvcheck 0.1.8 2020-03-01 [1]
S4Vectors * 0.24.4 2020-04-09 [1]
scales 1.1.1 2020-05-11 [1]
sessioninfo 1.1.1 2018-11-05 [1]
SRAdb 1.48.2 2019-12-24 [1]
stringi 1.4.6 2020-02-17 [1]
stringr 1.4.0 2019-02-10 [1]
SummarizedExperiment * 1.16.1 2019-12-19 [1]
survival 3.2-3 2020-06-13 [1]
sva * 3.34.0 2019-10-29 [1]
testthat 2.3.2 2020-03-02 [1]
tibble 3.0.3 2020-07-10 [1]
tidyr 1.1.0 2020-05-20 [1]
tidyselect 1.1.0 2020-05-11 [1]
usethis 1.6.1 2020-04-29 [1]
vctrs 0.3.2 2020-07-15 [1]
viridis * 0.5.1 2018-03-29 [1]
viridisLite * 0.3.0 2018-02-01 [1]
withr 2.2.0 2020-04-20 [1]
writexl 1.3 2020-05-05 [1]
xfun 0.16 2020-07-24 [1]
XML 3.99-0.3 2020-01-20 [1]
xml2 1.3.2 2020-04-23 [1]
xtable 1.8-4 2019-04-21 [1]
XVector 0.26.0 2019-10-29 [1]
yaml 2.2.1 2020-02-01 [1]
zlibbioc 1.32.0 2019-10-29 [1]
source
Bioconductor
Bioconductor
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
Bioconductor
Bioconductor
CRAN (R 3.6.1)
Bioconductor
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
Bioconductor
CRAN (R 3.6.1)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
Bioconductor
Bioconductor
CRAN (R 3.6.1)
Bioconductor
Bioconductor
Bioconductor
Bioconductor
Github (sinhrks/ggfortify@3f4020d)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
Bioconductor
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
Github (r-lib/memoise@4aefd9f)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
local
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
Github (easystats/report@dcdd283)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
Github (juba/rmdformats@94cd7a3)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.1)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.1)
Bioconductor
CRAN (R 3.6.2)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
Bioconductor
CRAN (R 3.6.2)
Bioconductor
[1] /home/ubuntu/R/x86_64-pc-linux-gnu-library/3.6
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library